Nonequilibrium steady state for strongly-correlated many-body systems: 

variational cluster approach 
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A numerical approach is presented that allows to compute nonequilibrium steady state properties 
of strongly correlated quantum many-body systems. The method is imbedded in the Keldysh 
Green's function formalism and is based upon the idea of the variational cluster approach as far as 
the treatment of strong correlations is concerned. It appears that the variational aspect is crucial 
as it allows for a suitable optimization of a "reference" system to the nonequilibrium target state. 
The approach is neither perturbative in the many-body interaction nor in the field, that drives the 
system out of equilibrium, and it allows to study strong perturbations and nonlinear responses of 
systems in which also the correlated region is spatially extended. We apply the presented approach to 
non-linear transport across a strongly correlated quantum wire described by the fermionic Hubbard 
model. We illustrate how the method bridges to cluster dynamical mean-field theory upon coupling 
two baths containing and increasing number of uncorrelated sites. 

PACS numbers: 71.27. +a 47.70.Nd 73.40.-c 05.60.Gg 



I. INTRODUCTION 



The theoretical understanding of the nonequilibrium 
behavior of strongly correlated quantum many-body sys- 
tems is a long standing challenge, which has become in- 
creasingly relevant with the progress made in the fields 
of quantum optics and quantum simulation, semicon- 
ductor, quantum, and magnetic heterostructures, nan- 
otechnology, or spintronics. In the field of quantum op- 
tics and quantum simulation recent advances in exper- 
iments with ultracold gases in optical lattices shed new 
light on strongly-correlated many body systems and their 
nonequilibrium properties. In these experiments, specific 
lattice Hamiltonians can be engineered and studied with 
a remarkable high level of control, making strong cor- 
relations observable on a macroscopic scale.— ~— In this 
field another very promising experimental setup to study 
correlation effects are coupled cavity quantum electrody- 
namic systems which contain some form of optical nonlin- 
earity resulting from the interaction of light with atomic 
levels 4^ These coupled cavity systems are inherently out 
of equilibrium, since they are driven by external lasers 
and susceptible to dissipation. Semiconductor, quantum, 
and magnetic heterostructures subject to a bias voltage 
also display nonequilibrium physics, where strong cor- 
relations play a decisive role. Experiments which study 
transport in molecular junctions demonstrate that many- 
body effects, also in combination with vibrational modes 
are crucial, see, e.g., Refs. 00. Another class of mate- 
rial structures with remarkable nonequilibrium proper- 
ties are (multi-well) heterostructures of diluted magnetic 
semiconductors (DMSs) and superlattices embedded in 
normal metals. These systems are of great interest as 
they open the possibility to tailor electronic and spin- 
tronic devices for computing and communications based 
on their unique interplay of spin and electronic degrees 
of freedom. Moreover, they display a pronounced non- 



linear transport behavior The source of nonlinearity 
is also related to the strong interaction between charge 
carriers, excitations and vibrational modes. In addition, 
spin degrees of freedom clearly play a major role in their 
transport properties. In order to fabricate technologi- 
cally useful structures the theoretical understanding of 
these highly correlated quantum many-body systems is 
indispensable. 

A typical nonequilibrium situation in all these sys- 
tems is conveniently described theoretically by switching 
on a perturbation at a certain time t = to, for exam- 
ple, a bias voltage, which is then kept constant after a 
short switching time. For this problem one may, on the 
one hand be interested in transient properties at short 
times after switching on the perturbation, for example 
in ultrafast pump-probe spectroscopy— In this case, the 
properties of the system depend on the initial state, as 
well as on the line shape of the switch-on pulse. For 
longer times away from To, quite generally one expects 
the system to reach a steady state, whose properties do 
not depend on details of the initial state. Nonequilib- 
rium steady states are relevant, for example, in quan- 
tum electronic transport across heterostructures, quan- 
tum dots, molecules (see, e.g., Refs. Il6l - [2ll) or in driven- 
dissipative ultracold atomic systemsi 22 Intriguingly, it 
was shown in Ref. [HI that nonequilibrium noise, which is 
present for instance in Josephson junctions, trapped ul- 
tracold polar molecules or trapped ions, still preserves 
the critical nonequilibrium steady states thus being a 
marginal perturbation as opposed to the temperature. 
Among the methods to treat strongly correlated sys- 
tems out of equilibrium, one should mention density- 
matrix renormalization group and related matrix- 
product state methods j 2 ^— continuum-time quantum 
Monte-Carlo,— different numerical and semi-analytical 
renormalization-group approaches 1 2 1 1 35 1 36 equation-of- 
motion methods) 16 ' 19 , dynamical mean-field thcory^2r— 
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scattering Bethe AnsatZ ) 41 ' 42 and the dual-fcrmion 
approach^ Recently, Balzer and Potthofi 4 ^ have pre- 
sented a generalization of cluster-perturbation theory 
(CPT) to the Keldysh contour, which allows for the 
treatment of time-dependent phenomena. Their results 
show that CPT describes quite accurately the short and 
medium-time dynamics of a Hubbard chain. A detailed 
study of the short-time dynamics of weakly correlated 
electrons in quantum transport based on the time evo- 
lution of the nonequilibrium Kadanoff-Baym equations, 
where correlations are treated in Hartree-Fock-, second 
Born-, and GW- approximation has been given in Ref . l45l. 
These approximations are restricted to moderate cor- 
relations but on the other hand they allow to study 
rather complex models and geometries. As far as the 
steady-state behavior is concerned, the nonequilibrium 
(Keldysh) Green's function approach has been widely 
used on an ab-initio or tight-binding level, where correla- 
tions are treated in mean-field approximation. Since the 
effective particles are non-interacting, the Meir-Wingreen 
expression 1 ^ for the current can be applied, which re- 
lates the current to the retarded Green's functions of the 
scattering with a self-energy that is renormalized due to 
the presence of the leads. Representative applications 
for nano-structured materials and molecular devices are 
given in Refs. I46U481 and in the review article Ref. [2(| 

Here we aim at strongly correlated many-body sys- 
tems, and we propose a variational cluster method, that 
allows to study steady-state properties. 

The paper is organized as follows: In Sec. [H] we present 
the variational cluster method to treat correlated systems 
out of equilibrium. After an introductory discussion as 
well as relation to previous work, we present the general 
method in Sec. Ill Al We discuss the self-consistency con- 
dition in Sec. Ill Bl In Sec. Mil we introduce two specific 
models describing a strongly correlated Hubbard chain 
and a strongly correlated Hubbard ladder, respectively, 
which are embedded between left and right uncorrelated 
reservoirs with different chemical potentials and on-site 
energies. This results in a voltage bias which is applied to 
the system. Results for the steady-state current density 
are discussed in Sec. IIVI Finally, in Sec. [V] we present 
our conclusions and outlook. 



II. METHOD 

In order to study nonequilibrium properties of strongly 
correlated systems one typically considers a model con- 
sisting of two leads with uncorrelated particles, and a 
central correlated region. The three regions are initially 
decoupled. At a certain time tq a coupling V between the 
three regions is switched on. A natural approach is to 
treat V via strong-coupling perturbation theory, which 
at the lowest order essentially corresponds to cluster- 
perturbation theory (CPT) . In Ref. |3 it has been shown 
that the short time behavior can be well described within 
CPT. This can be understood from the observation that 



switching on the inter-cluster hopping V for a certain 
time At produces a perturbation of order V At, which 
is accounted for at first order in CPT. Therefore, we 
expect the result to be accurate for small At. When 
addressing the steady state it is, thus, essential to im- 
prove the long-time behavior. Here, we suggest that 
nonequilibrium CPT can be systematically improved by 
minimizing some suitable "difference" between the un- 
perturbed ( "reference" ) state which enters CPT and the 
target steady state. 

The strategy presented here to achieve this goal con- 
sists in exploiting the fact that the decomposition of the 
Hamiltonian into an "unperturbed part" and a "per- 
turbation" is not unique. Prompted by the variational 
cluster approach (VCA), one can actually add "auxil- 
iary" single particle terms to the unperturbed Hamilto- 
nian and subtract them again within CPT. This freedom 
can be exploited in order to "optimize" the results of 
the perturbative calculation. As discussed in detail in 
Refs. I49ll50l in equilibrium this is an alternative way to 
motivate the introduction of variational parameters in 
VCA. The idea discussed here, thus, provides the natu- 
ral extension of VCA to treat a nonequilibrium steady 
state. There remains to define a criterion for the "dif- 
ference" between initial and final state. (Cluster) Dy- 
namical Mean-Field Thcory^^ 1 .^ (DMFT) provides a 
natural solution, requiring the cluster-projected Green's 
functions of the initial and final state to coincide. Of 
course, this self-consistency condition requires an infinite 
number of variational parameters, as well as the solution 
of a (cluster) impurity problem, which is computationally 
very expensive and whose accuracy is limited, especially 
in real time. In equilibrium, the self-energy functional 
approac h 54 ! 55 (SFA) provides one possible generalization 
of DMFT if one wants to restrict to a finite number of 
variational parameters. In this case, the requirement for 
the "difference" is provided by the Euler equation (see, 
e.g., Eq. (7) in Ref-H. 

In the present paper, we explore an alternative cri- 
terion, represented by (fT5)l . which, upon including an 
infinite number of bath sites, becomes equivalent to 
(cluster)-DMFT (see App.E}, similarly to SFAM With- 
out bath sites this corresponds to requiring that, for a 
given set of variational parameters p, their conjugate op- 
erators, i.e., dh/dp, h being the Hamiltonian, have the 
same expectation value in the unperturbed and in the 
final target state. This criterion is numerically easier to 
implement than the SFA, since in this case it is not neces- 
sary to search for a saddle point, which is well known to 
be numerically expensive^ In addition, inclusion of bath 
sites provides self consistency conditions for dynamic cor- 
relation functions as well. 

The freedom discussed above can be additionally ex- 
ploited by including the hybridization between correlated 
regions and the leads as well as part of the leads them- 
selves into the unperturbed Hamiltonian which is solved 
exactly by Lanczos exact diagonalization. In this way, 
CPT is then used to treat hopping terms further away 
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from the correlated region^ This partly accounts for the 
influence of the leads onto the self-energy of the corre- 
lated region. 

Finally, let us mention that the method is probably 
most suited to deal with models for which the correlated 
region is spatially extended (see Fig. [I}. In this case, 
this region must be partitioned into clusters which can be 
solved exactly, while the intercluster terms are included 
into the perturbativc part. 

A. Variational cluster approach for nonequilibrium 
steady state 

The physical model of interest consists of a "left" and 
"right" noninteracting lead, as well as a correlated region 
described by the Hamiltonians hi, h ri and h c , respec- 
tively, see Fig.[T] h c contains local (Hubbard-type) in- 
teractions, as well as arbitrary single-particle terms. For 
r < To, the three regions are in equilibrium with three 
reservoirs at different chemical potentials, fii, fi r , and fi c 
respectively. The correlated region is much smaller in size 
than the leads, so that the latter act as relaxation baths. 
At t = To, the single particle (i.e., hopping) Hamilto- 
nian terms Vi c and V rc are switched on. These connect 
the left and right reservoir, respectively, with the corre- 
lated region. The total time-dependent Hamiltonian is, 
thus, given by 

H(T) = h + 6(r-T )T , (1) 

where h = h c + hi + h ri and T = Vi c + V rc . We con- 
sider here the fermionic case, although many concepts 
can be easily extended to bosons. After a time At long 
enough for relaxation to take place, the system reaches a 
nonequilibrium steady-state, with a particle current flow- 
ing from left to right for fii > fi r and from right to left 
for hi < \i r . 

As discussed above, the total r > tq Hamiltonian H = 
H(t > To) is decomposed into an unperturbed part h 
and a perturbation T: 

H = h + f. (2) 

In the simplest CPT approach for a "small" correlated 
region one can take h = h, and T = T. However, when 
the correlated region is extended, as in Fig. [TJ it has 
to be further decomposed into smaller clusters that can 
be solved by exact diagonalization^ In this case, the 
intercluster hopping is subtracted from h and must be 
included in T. In addition, one can include part of the 
leads into the clusters (dashed lines in Fig. [1]), so that 
Vic + V rc are incorporated into h, while the leads inter- 
cluster hoppings (e.g. tu c in the figure) are included-^ 
in T. Finally, in the spirit of VCA, arbitrary intracluster 
terms Ah can be added to the unperturbed Hamiltonian 
and subtracted perturbativcly within T. In other words, 
calling h c i the Hamiltonian describing the physical clus- 
ter partition, and T c i the one describing the intercluster 
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FIG. 1: Generic scheme of the model studied here: full 
(empty) circles indicate correlated (uncorrelated) lattice sites. 
Correlated sites define the correlated region (c), and are char- 
acterized by an on-site Hubbard interaction U, an on-site en- 
ergy e c , and by hopping elements t x and t y in the x and y 
direction, respectively. The physical leads (l,r), indicated by 
the two shaded areas, consist of half-infinite planes described 
by uncorrelated tight-binding models with hopping tL, on-site 
energies ei and e r , and chemical potentials fii and /i r , respec- 
tively. The correlated region is connected to the leads via hop- 
pings V. The width (number of sites in the x direction) of the 
correlated region is L cx . The height of the whole system in 
the y direction is infinite. In this work, we study two 
strongly correlated chain (L cx = 1) and a strongly correlated 
two-leg ladder (L cx = 2), both perpendicular to the applied 
bias. In the variational cluster calculation the central region 
described by the unperturbed Hamiltonian h can differ from 
the physical one. The latter coincides with the correlated sites 
(white area in the figure).— On the other hand, the former 
consists of disconnected clusters aligned along the y direction, 
one of them being represented by the dash-dotted rectangle 
in the figure. The corresponding equilibrium Green's func- 
tion is determined by Lanczos exact diagonalization. The 
size of these clusters is L c = L cx x L cy (4 x 2 in the example) . 
The coordinates of the left and right boundary sites of the 
central region are indicated by xu and x^ T , respectively. Ac- 
cordingly, dashed lines represent hopping processes, which are 
omitted in the unperturbed (reference) Hamiltonian h and are 
re-included perturbatively within T. Full lines indicate hop- 
ping terms present in h, which are thus treated exactly (see 
text). 

hoppings (dashed lines in Fig. [T]), we write h = h c i + A/;, 
and T = T c i — Ah so that the total Hamiltonian remains 
unchanged: 

H = hd+Tci = h + f . (3) 

The arbitrariness in the choice of Ah can be exploited to 
optimize the unperturbed state, as discussed in Ref. [5^ 
for the equilibrium case. Here, we will adopt a differ- 
ent optimization criterion, see discussion below. Being a 
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single-particle term, T is described by its hopping ma- 
trix T. This matrix has a block structure according to 
the three regions discussed above and shall be denoted 
by Tic, T rc and T cc , respectively. 

Nonequilibrium properties, in general, and nonlinear 
transport in particular can quite generally be deter- 
mined in the frame of the Keldysh Green's function ap- 
proachi 16 i 17 i 59 ~ — Here, we adopt the notation of Ref. Il7l . 
for which the 2x2 Keldysh Green's function matrix is 
expressed as 

G(v,v'\t,t')=(C R ^) , (4) 

where the retarded (G R ), advanced (G A ), and Keldysh 
(G K ) Green's functions depend in general on two lat- 
tice sites (r, r') and two times (r, r'). However, both for 
t < To as well as in steady state, time translation invari- 
ance holds, so that Green's functions depend only on the 
time difference t — t', and we can Fourier transform to 
frequency space uj. 

We use uppercase letters G to denote Green's functions 
of the full Hamiltonian H, and lowercase g for the ones 
of the unperturbed Hamiltonian h. The advantage of 
using the Keldysh Green's function matrix representation 
is that one can express Dyson's equation in the same form 
as in equilibrium! 16 ' 17 In our case, we can express it in 
the form 

G = g + g (T + AS) G , (5) 

where g = diag(<?;z, gcc, grr) is block diagonal, and the 
products have to be considered as matrix multiplica- 
tions^ In ([5]), AS = S — S/j is the difference between 
the (unknown) self-energy S of the total Hamiltonian H, 
including the coupling to the leads, and the self-energy 
Hh associated with the unperturbed Hamiltonian h. 

The CPT approximation 6 ^ precisely amounts to ne- 
glecting AS. As pointed out in Ref. HfH this corresponds 
to neglecting irreducible diagrams containing interactions 
and one or more T terms. It should, however, be stressed 
that the self-energy of the isolated clusters is exactly in- 
cluded in g cc , which is obtained by Lanczos exact diago- 
nalization. 

In this approximation, ([5]) can be used to obtain an 
equation for the Green's function G cc projected onto the 
central region, which is still a matrix in the lattice sites 
of the central region and in Keldysh spaced (this is a 
straightforward generalization of, e.g., the treatment in 
Ref. til): 

e{Lr} 

Gcc = 9cc + 5cc( T cc G cc + ^ ' T ca G a c) (6) 

a 

and for the lead-central region Green's functions: 

Gac — g<xa T ac G cc , with a 6 {l,r}. (7) 



It is noteworthy that Eq. ([7| is exact and not based 
on the CPT approximation, as the leads contain non- 
interacting particles. Insertion of (O into ^ yields 

G C c = gcc + 9cc\T C c + S cc ) G cc (8) 
with the lead-induced self-energy renormalization 

E{l,r} 

Sec — ^ T C ct gaa T a c ■ (9) 
a 

Here g aa stands for the Green's function of the isolated 
lead a. One finally obtains a Dyson form for the steady 
state Green's function of the coupled system at the cen- 
tral region 

Gcc = 9cc - Tcc - ^ ■ (10) 

Different from the usual Dyson equation, g cc is the 
Green's function for the isolated clusters, which contains 
all many-body effects inside the cluster. 

For the evaluation of the current from, say, the left 
lead to the central region^ 7 - one needs the G; c Green's 
function, which is readily obtained by combining ([7]) 
with ([TO]) . This leads to the generalized Kadanoff-Baym 
equation (see e.g. Refs. [l6lll8h . along with the fact that 
the central region is finite in x direction and the leads 
are infinite, one can rewrite the current into a Biittiker- 
Landauer type of formula 

i = J ~ Mr) - fp(e - /i;)] 

xTr[ G£(e)r,( e ) Gf c (s)T r (s)] . (11) 

where G^J A is the retarded/advanced part of the Green's 
function G cc , and the trace, as well as matrix products 
run over site indices in c. T a describes the inelastic 
broadening owing to the coupling to lead a, which in 
CPT is given by 6 ^ 

T Q = 21m {T ca g aa Tac} , 

which represents the contribution of lead a to the imag- 
inary part of S^,. Interestingly, the expression for the 
current in CPT has the same structure as the Meir- 
Wingreen formula— for non-interacting particles, which 
is the basis for nonequilibrium ab-initio-calculations, 47 
Here, however, the Green's function contains the many- 
body interactions of the correlated region. An advantage 
of this expression is that it yields an explicit connection to 
the Green function G^I A of the scattering region and the 
influence of the itinerant electrons in the leads. A simi- 
lar expression can be derived for the one-particle density 
matrix between two sites with the same y coordinate, 
which is required for the self-consistency condition dis- 
cussed below. 

As it is well known, all retarded and advanced Green's 
functions are evaluated without chemical potentials. The 
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latter enter through the Keldysh Green's function or 
rather via the Fermi functions. While the chemical po- 
tential of the central region is wiped out in the steady 
state due to its small size in comparison to the size of the 
leads, the chemical potentials of the leads explicitly enter 
the expressions for the current and the density matrix, 
see Eq. (fTTT) . In the case investigated here, the central 
region is translation invariant in y direction and is split 
into identical clusters. In the end, as far as the main 
numerical task is concerned, one has to solve many-body 
problems for clusters of size L = L cx x L cy , invert ma- 
trices of the same size, and sum over wave vectors q y 
belonging to the Brillouin zone associated with the clus- 
ter supercell. 



to noncquilibrium cluster DMFT. Generalization of the 
SFA condition to nonequilibrium should be, in principle, 
obtained by replacing gocc with £/, in (fT3|) . 

A second systematic improvement of this nonequilib- 
rium VCA approach consists in increasing the cluster size 
L c . This can be done in two ways: (i) by extending the 
boundaries of the central region in y direction and thus 
treating more correlated sites exactly^ and (ii) by ex- 
tending the boundaries in x direction to include an in- 
creasing number of uncorrelated lattice sites, i.e., taking 
L cx > L CXl cf. Fig. [T] This amounts to taking into ac- 
count to some degree the ^-induced renormalization of 
the self-energy. 



B. Self-consistency condition 



Equation (|10[) is the expression for the Green's func- 
tion of the central region within the CPT approximation. 
As discussed above, one would like to optimize the ini- 
tial state in some appropriate way by suitably adjust- 
ing the parameters Ah of the unperturbed Hamiltonian 
h. The inclusion of additional terms Ah adds flexibil- 
ity to the self-energy E/j which is included within this 
approximation. Obviously, it makes no difference in the 
case of non-interacting particles as the selfenergy van- 
ishes exactly, independently of Ah. This freedom can be 
exploited in order to improve the approximation system- 
atically. A similar discussion on this issue has been given 
in Refs. SH3), and is at the basis of the VCA ideaM. 

As discussed above, we need a variational condition as- 
sociated with a "minimization" of the difference between 
unperturbed and perturbed state. In (cluster)-DMFT 
one requires the cluster projected Green's function to be 
equal to the unperturbed one 



g cc = V(G CC ) 



(12) 



where V projects the Green's function onto the cluster, 
i. e., it sets all its intercluster matrix elements to zero.— 
Since here we have a finite number of variational param- 
eters p that can be adjusted, we cannot satisfy (|T2|) . We, 
thus, propose a "weaker" condition, namely that the ex- 
pectation values of operators coupled to the variational 
parameters contained in Ah (i.e., dAh/dp) be equal in 
the unperturbed and in the perturbed state. More specif- 
ically, we impose the condition 



du; „ d(g 0cc )~ 
2^ trTl ^p- 



(g cc - G cc ) = , 



(13) 



where t\ is a Pauli matrix in Keldysh space,— and gocc 
is the Green's function associated with the noninteract- 
ing part of h. It is interesting to note (see appendix 
|A")) that by including into Ah a coupling to an infinite 
number of bath sites, the present method, with the self- 
consistence condition ()13|) whereby p are the bath param- 
eters (hopping and on-site energies), becomes equivalent 



III. MODEL 

Next, we present an application of the nonequilibrium 
VCA method described in Sec. HU Specifically, we study 
nonlinear transport properties across an extended corre- 
lated region (denoted as c in FigJlJ, which we take to be a 
Hubbard chain (L cx = 1) or a Hubbard ladder (L cx = 2) 
with nearest-neighbor hoppings t x and t y , on-site inter- 
action U, on-site energy e c , and chemical potential fi c 



tij C i<j C ja 



U ^2 n^hix + (e c - /i c ) ^ 



in usual notation, and where ty = t x (ty = t y ) for i and 
j being nearest neighbors in x direction (y direction). 
The leads (shaded regions in Fig. [T]) are described by 
two-dimensional semi-infinite tight-binding models with 
nearest-neighbor hopping t^, on-site energies e/ and e r , 
and chemical potentials fii and /i r for the left and right 
lead, respectively. We apply a bias voltage V& to the leads 
by setting \i T = e r = V&/2 and restrict to the particle-hole 
symmetric case where e c = —U/2, (x c = 0, e r = — e;, and 
111 = —fir. For simplicity, we neglect the long-range part 
of the Coulomb interaction. Under some conditions, this 
can be absorbed within the single-particle parameters of 
the Hamiltonian, in a mean-field sensed 

As discussed above, the unperturbed Hamiltonian h 
does not necessarily coincide with the physical partition 
into leads and correlated region, h is obtained by tiling 
the total system into small clusters as illustrated in Fig.[TJ 
as well as by adding an intracluster variational term Ah. 

In the present work Ah describes a correction At x 
to the intra-ladder hopping. Further options could in- 
clude, for instance, a site-dependent change in the on- 
site energy Ae c (x). Particle- hole symmetry can be pre- 
served by constraining this change to be antisymmetric: 
Ae c (x) = — Ae c (— x). In this paper, whose goal is to 
carry out a first test of the method, we restrict, for sim- 
plicity, to a single variational parameter. The choice of 
At x as a variational parameter is motivated by the fact 
that this term is important for the current flowing in x 
direction. According to the prescription discussed above, 
we require the expectation value of the one-particle den- 
sity matrix for nearest-neighbor indices in x direction to 
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FIG. 2: (Color online) Steady-state current density j x versus bias voltage VJ, for a correlated two-leg ladder (L cx — 2). First row 
shows j x normalized by V 2 as function of VJ, evaluated for different values of V and of the interaction (a) U = 2.0, (b) U = 4.0, 
and (c) U = 6.0. Second row shows the U dependence of the current for different values of the hopping V = Vi c = V TC from the 
leads to the correlated region (d) V = 1.0, (e) V = 0.5, and (f) V = 0.1. Solid (dashed) lines represent results for the current 
between the left lead and the central region (between two in x direction adjacent sites inside the central region), i. e., evaluated 
with (G^), see text for details. Results are obtained by using a reference Hamiltonian h consisting of disconnected clusters 
of size L c = L cx x L cy = 2x6. 



be the same for the unperturbed h and for the full H, 
i.e. evaluated with g cc and with G cc . 

One comment about the chemical potential. In princi- 
ple, when including some of the sites of the leads in h, i. e., 
when L cx > L cx , then these additional sites have a chem- 
ical potential fi c which differs from the one they would 
have if L cx = L cx (i.e., /i; or fx r ). However, the chemical 
potential, of these sites does not affect the steady state, 
as their volumc-to-surface ratio is finite. Of course, their 
on-site energies (e r and ej) are important. 

Due to translation invariancc by a cluster length L cy in 
the y direction, it is convenient, as in usual VCA, to carry 
out a Fourier transformation in y direction, with associ- 
ated momenta q y . The Green's functions g cc and G cc , as 
well as T become now functions of two momenta q y + Q y 
and q y + Q' yl where Q y and Q' y are reciprocal super lattice 
vectors of which there are only L cy incquivalent ones. In 
order to evaluate the nonequilibrium steady state, one 
only needs the equilibrium Green's function g(xb a \qy\z) 
of the isolated leads at the contact edge to the central 
region, with x coordinate equal to Xb a (a € {/, r}), and 
Fourier transformed in the y directions, where q y is the 
corresponding momentum and z the complex frequency 
For a semiinfinite nearest-neighbor tight-binding plane 
with hopping t^, and on-site energy e Q , this can be ex- 



pressed as 

g{x ba \q y \z) = g c ,ioc{z - 2t L cosq y - e a ) , (14) 

where g c .ioc(z) is the local Green's function of a tight 
binding chain with open boundary conditions and with 
zero on-site energy. The latter can be determined ana- 
lytically along the lines discussed in Ref. |68|. 

The model studied here, is motivated by the interest 
in transport across semiconductor heterostructures (see, 

e.g. miisizjj). However, it is well known that in this 

case charging effects arc important, also near the bound- 
aries between the leads and the central region. Here, scat- 
tering effects produce charge density waves, which, when 
taking into account the long-range part of the Coulomb 
interaction, even in mean-field, produce a modification 
of the single-particle potential. In order to treat real- 
istic structures, these effects should be included at the 
Hartree-Fock level at least. All these generalizations can 
be straightforwardly treated with the presented varia- 
tional cluster method, however, in this work we focus on 
a first proof of concept study and application contain- 
ing the essential ingredients for the investigation of the 
nonequilibrium steady state of strongly correlated many- 
body systems. 
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FIG. 3: (Color online) Steady-state current density j x as in Fig. [2] but for the correlated chain {L cx — 1). The current density 
is evaluated for different values of the lead to correlated region hopping (a) V = 1.0, (b) V = 0.5, and (c) V — 0.1, and of the 
interaction U, see legend. Results are obtained for reference clusters of size L c = L cx x L cy = 3x4. 



IV. RESULTS 

We have evaluated the steady-state current density j x 
of the models discussed in Sec. Mil as a function of the 
bias Vb = e r — e; between the leads at zero temperature. 
Simultaneously the chemical potential is adjusted to the 
on-site energy /i Q = e Q , which corresponds to a rigid 
shift of the density of states in both leads in opposite 
directions. 

In Fig. [5] we display results for the two-leg ladder 
(L cx = 2), for different values of the interaction strength 
U = {0,2,4,6} and lead-to-system hopping V = 
{1.0, 0.75, 0.5, 0.25, 0.1}. We use h = 1 and t L = 1 
which sets the unity of energy. Moreover, we take the lat- 
tice constant a = 1 . The hopping is uniform in the whole 
system, meaning that t x , t y in the correlated region and 
iz, of the leads are equal. The on-site energy of the cor- 
related region is e c = — U/2 corresponding to half- filling, 
whereas the on-site energy of the left (right) lead is equal 
to its chemical potential fii (fi r ). The unperturbed hamil- 
tonian h describes the central region decomposed into 
clusters of size L c = 2 x 6. The corresponding Green's 
function g cc is determined exactly by Lanczos diagonal- 
ization. All results are determined self-consistently using 
At x as variational parameter, see Sec. IIIII 

Using the Meir-Wingreen expression, Eq. (fTTj) . the gen- 
eral trend of the results for the steady-state current j x 
can be discussed conveniently. At zero temperature there 
are only contributions to the current for min(//z,/v) < 
u < max(/ii, /z r ) due to the difference of the Fermi dis- 
tribution functions. In particular this leads as expected 
to zero current for zero bias voltage Vb- With increasing 
bias voltage Vb the modulus of j x initially increases. For 
large values of Vb it decreases again, as the overlap of the 
local density of states of the two leads enters the expres- 
sion, which is zero if Vb is greater than the band width of 
the leads. Hence the local density of states of the leads 
along with the Fermi function act as a filter that averages 
the electronic excitations of the central region within a 
certain energy window. 



In the system we are studying, the leads are modeled 
by semi-infinite tight binding planes. Alternatively, in- 
stead of using (fl~4"|) one could simply put a model Green's 
function "by hand," as for example one which describes 
a Lorentzian shaped density of states. Such an unbound 
density of states generally leads to a finite value of the 
current for arbitrary bias. 

The leads have a further effect on the result as they 
provide an inelastic broadening of the energy spectrum 
of the central region entering E eff , see Eq. (fl"0)). which 
smears out details of the excitation spectrum. As far 
as the lead-correlated region coupling V is concerned, 
there are two competing effects: on the one hand, the 
current increases with increasing V due to the stronger 
coupling between the correlated region and the leads. On 
the other hand, details of the electronic excitations are 
smeared out with increasing V leading to a reduced reso- 
lution. Therefore, in order to detect the effects of strong 
correlations, particularly the gap, a small value for V is 
required. 

The details of the V dependence of j x for small V 
can be deduced from (fTTj) . Here, the expression for the 
current has a prefactor proportional to V A (at least in 
the L cx = L cx case), due to the two T terms. On the 
other hand, for a gapless system, there is a V 2 term in 
the denominator of \G^ C \ 2 . For a gapped system, this 
is cut off by the energy gap E g , so that in this case 
jx ~ V^/Eg, while j x ~ V 2 for a gapless spectrum. 
These aspects are clearly observable in Fig. [5] (a)-(c), 
which shows the scaled current density j x /V 2 for fixed 
interaction strength U but varying V. The envelope has 
a rotated S-like structrue due to the combined effects of 
the lead density of states and of the Fermi functions. 

Next we will analyze a bit more in detail the effects 
of the Hubbard interaction. Increasing the interaction 
strength U in the correlated region leads to a suppression 
of the current and the opening of a gap, which is best 
oberserved in (f) . For U = 4 the maximum of the current 
density is roughly reduced by a factor of two as compared 
to the noninteracting case, whereas for U = 6 the current 
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FIG. 4: (Color online) Convergence 
of the steady-state current density 
j x with reference cluster size L c — 
L cx x L cy for the correlated two-leg 
ladder with V = 0.5. Results in 
(a), (b) are obtained by a variational 
adjustment of the intra-cluster hop- 
ping t x as discussed in the text, while 
those of (c), (d) are obtained without 
modification of t x . The values for 
the Hubbard interaction are U = 2 
in [(a), (c)] and U = 4 in [(b), (d)]. 



is almost one order of magnitude smaller as compared to 
the noninteracting system, see Fig. [2] (d)-(f). 

Finally, we want to address the the difference between 
the solid lines and dashed lines in the the panels (d)-(f) 
of Fig. [21 which represent the current density evaluated 
on a bond connecting the leads to the central region, or 
on a bond within the two-leg ladder. Due to the sta- 
tionary condition, the two results should coincide. How- 
ever, our calculations shows a slight discrepancy between 
them, which is due to the fact that the method is not 
completely conserving and, thus, the continuity equation 
is not completely fulfilled. However, from our results we 
see that the deviation from the continuity equation is 
quite small. We expect this discrepancy to be reduced 
upon improving the optimization with the introduction 
of additional variational parameters. 

In Fig. [3] we show the steady-state current density j x 
across the correlated chain {L cx = 1) as a function of 
the bias voltage. The parameters are the same as in the 
case of the two-leg ladder, however, the central region 
is decomposed into clusters of size L c = 3x4, where 
also sites of the leads are taken into account to improve 
the results. The half-filled Hubbard chain is gapped as 
well. As for the two-leg ladder, the gap behavior can 
be better seen in the current-voltage characteristics for 
smaller values of V, in our case for V = 0.1. In contrast, 
for strong coupling V = 1.0, (a), no gap behavior can be 
seen in the current due to the strong hybridization with 
the leads. 

For strong values of the coupling V between leads and 
correlated region (V = 1.0), (a), the current is signifi- 
cant for all values of the interaction U. However, with 
decreasing V , (b)-(c), the current is strongly suppressed 



for large interaction U. Importantly, for the correlated 
chain the continuity equation is always strictly fulfilled. 
In other words, there is no difference between j x evalu- 
ated on a intercluster bond between the leads and the 
cluster, or on an intracluster bond. This is due to the 
absence of vertex corrections at the uncorrclated sites. 

Next, we study the convergence of our results with 
the size of the cluster, as well as the effect of the self- 
consistency condition for the two- leg ladder and V = 0.5. 
Results are depicted in Fig. @] for two different values of 
the Hubbard interaction, namely U = 2 [(a), (c)] and 
U = 4 [(b), (d)]. We do not plot results of the con- 
vergence analysis for U — 6, since for this large U the 
current is already rather small, as can be seen in Fig. [2] 
(d)-(f) . Results in (a) and (b) , first row, are obtained by 
adjusting At x self-consistently, as described in Sec. IIII1 
whereas (c) and (d), second row, shows results without 
self-consistency, i.e., with At x = 0. Results show that 
the self-consistency procedure improves the results, as 
the convergence for j x is faster with increasing cluster size 
as compared to the case without self-consistency. Gen- 
erally, we observe pronounced finite size effects for very 
small clusters up to 2 x 4, and convergence seems to be 
reached for the 2x6 cluster. 

We now repeat the same analysis for the correlated 
chain. The corresponding current densities for the pa- 
rameters U — 2 and V = 0.5 are shown in Fig. [5] for 
different cluster sizes. Results shown in (a) are with self- 
consistency procedure (flU)) . whereas the results shown 
in (b) are without. In the present case, where we con- 
sider transport across a strongly correlated chain, con- 
vergence is achieved very quickly with increasing cluster 
size. Therefore, there is no sensible difference between 
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FIG. 5: (Color online) Convergence of the steady-state cur- 
rent density j x with reference cluster size L c — L cx x L cy for 
the correlated chain. Results in (a) fulfill the self-consistency 
condition (| 13f) . whereas results in (b) do not. The parameters 
are U = 2 and V = 0.5. 



results obtained with or without self-consistency, apart 
for the pathological case L c = 3 x 1 (see below). 

Results obtained for the two-leg ladder and for the 
chain show that cluster geometries with L cy = 1 provide 
results far from convergence, even with self-consistency. 
For the chain this is probably due to the degeneracy of 
the cluster ground state. For the ladder, it seems that 
using as starting point the 2x1 dimer exaggerates the 
gap. But besides these data obtained from admittedly 
very small clusters, results converge quickly MS (\ func- 
tion of cluster sizes, especially when the hopping in x 
direction is used as a variational parameter. 



V. CONCLUSIONS 

In this paper we have presented a novel approach to 
treat strongly correlated systems in the nonequilibrium 
steady state. The idea is based on the variational cluster 
approach extended to the Keldysh formalism. For the 
present approach the expression for the current resembles 
the corresponding Meir-Wingreen formulas. As in the 
original Meir-Wingreen approach, which is also the basis 



for nonequilibrium density-functional based calculations, 
we directly address the steady state behavior of a device 
coupled to infinite leads. The latter is necessary for the 
system to reach a well-defined steady state. 

The present nonequilibrium extension is in a similar 
spirit to the equilibrium self-energy functional approach, 
in which one "adds" single-particle terms to the cluster 
Hamiltonian which is then solved exactly, and "subtract" 
them perturbatively i 49 ' 50 The values of the parameters 
are determined by an appropiatc requirement which in 
the end amounts to optimizing the unperturbed state 
with respect to the perturbed one. 

There is a certain freedom in choosing the most appro- 
priate self-consistency criterion. Here we have required 
the operators associated with the variational parameters 
to have the same expectation values in the unperturbed 
and in perturbed state. Certainly, an interesting alter- 
native would be to generalize the variational criterion 
provided by the self-energy functional approach^ to the 
nonequilibrium case. This will be obtained by a suit- 
able generalization of the Euler equation (Eq. (7) of 
Ref. |54[ ) to the Keldysh contour, i.e., by replacing go cc 
with the self-energy £/, in (fT3"|) Work along these lines is 
in progress. 

The advantage of the present variational condition (|13p 
is that it is computationally less demanding, as one just 
needs to evaluate cluster single-particle Green's func- 
tions. Which one of the two conditions gives more ac- 
curate results cannot be stated a priori and should be 
explicitly checked. 

In any case, both methods, the self-energy functional 
approach and the present one, become equivalent to 
(cluster) dynamical mean-field theory in the case in 
which an infinite number of variational parameters is 
suitably taken (see Appendix [A}. 

In general, we expect results to improve when more 
variational parameters are taken into account. In partic- 
ular, when evaluating the current across the central re- 
gion, it would be useful if a current was already flowing 
in the cluster. This can be achieved by adding a complex 
variational hopping between the end points of the cluster, 
and of course remove it perturbatively. The correspond- 
ing variational condition would contain the interesting 
requirement that the current flow in this modified clus- 
ter be the same as in steady state. 

The model studied here, is motivated by the interest 
in tr anspor t across semiconductor heterostructures (see, 

e.g. yiijiiizO). However, it is well known that in 

this case charging effects are important, also near the 
boundaries between the leads and the correlated region. 
Here, scattering effects produce charge density waves, 
which, when taking into account the long-range part of 
the Coulomb interaction, even in mean-field, produce a 
modification of the single-particle potential. In order to 
treat realistic structures, these effects should be included 
at the Hartree-Fock level at least. All these general- 
izations can be straightforwardly treated with the pre- 
sented variational cluster method, however, in this work 
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we focus on a first proof of concept study and application 
containing the essential ingredients for the investigation 
of the noncquilibrium steady state of strongly correlated 
many-body systems. 
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Appendix A: Connection to (cluster) Dynamical 
Mean-Field Theory 

Here, we show that the self-consistent condition ([15)) 
provides a bridge to (cluster) DMFT, when an increas- 
ing number of noninteracting bath sites with appropri- 
ate parameters and occupations is included in the clus- 
ter Hamiltonian. Notice that these are "auxiliary" baths 
and are not related to the leads. Concretely, this is 
achieved by introducing into the variational Hamiltonian 
Ah a coupling of the central region with a set of uncorre- 
cted bath sites with appropriate energies, hybridizations 
v^, and occupations (see below). The hybridizations 
and the energies are therefore "included" in h, but "sub- 
tracted perturbatively," from the target Hamiltonian H. 
Their parameters are determined variationally via (|13|) . 

Now, since go cc is cluster-local, a solution to (|13[) is 
obviously given by (|12[) . However, this solution can gen- 
erally not be obtained with a finite number of parameters. 
As in usual equilibrium (cluster) DMFT^ (fT2")) can thus 
be solved via an iterative procedure defined by 

9occ,new = (P(Gcc))~ + 

= 9occ,old ~ 9cc ■ (Al) 

It is, thus, sufficient to show that an arbitrary gocc,new 
can be obtained by coupling the cluster to a noninter- 
acting bath with suitably chosen bath parameters. For 
the retarded and advanced Green's functions, the proce- 
dure is the same as in equilibrium. The Keldysh part 
is slightly more complicated. In order to show that an 
aribrary go cc ,new can be realized, one introduces the hy- 
bridization function 



where the A R , A A , and are matrices in the clus- 
ter sites. Similarly to equilibrium DMFT gocc,new is ex- 
pressed as 

90cc.new = So~cc,0 ~ A M ■ ( A 3) 

Here, g 0cc is the "bare" noninteracting cluster Green's 
function, i. e., the one with neither baths nor variational 
parameters. 

An arbitrary (steady-state) A(w) can be produced by 
coupling the central region to an appropriate bath in the 
following way. The retarded (and advanced) part are 
obtained as in equilibrium DMFT— by coupling to a bath 
with spectral function^ A bat h{oj) given by 

A bath {u) = --lmA R {u) , (A4) 

(ReA 71 is fixed by Kramers-Kronig relations). On the 
other hand, the Keldysh part is generated by splitting 
the bath defined by (|A4[) into two baths, a full (fj, — 
oo) and an empty (p = — oo) one, respectively. Their 
spectral functions are denoted by A^ ath (uj) and A^ ath {ui), 
respectively, and should obviously fullfill the condition 

A bath {u)=A+ ath + Ai ath . (A5) 

Since the Fermi functions of the two baths are 1 and 0, 
respectively, the Keldysh part A k (oj) is given by (A^ is 
anti-hcrmitian) 

A*(w) = -2iir B bath (to) ee -2in (A^ ath (cj) - A+ t/ ») . 

(A6) 

This fixes the two spectral functions to be 
, T , s _ Abath{u) ± B bath (u>) 

A batfiM - 2 ' ^ ' 

As usual, the two baths spectral functions A bat h,±(w) are 
realized by coupling the central region with a set of non- 
interacting sites with energies and hybridizations^ 
v^, fixed by 

At th {u) = Y, v t<^("-zt)- (A8) 

n 
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